Evaluación 3¶
Patiy Li Yang
Ejercicio 1: Frontreras de China¶
In [1]:
import os
import geopandas as gpd
from fiona import listlayers
ChinaMaps_3415="https://github.com/Ciencia-de-datos-espaciales-2023-2/Evaluacion3_GeoDF-Operaciones/raw/main/maps/China3415/ChinaMaps_3415.gpkg"
China3415 = gpd.read_file(ChinaMaps_3415,layer="countries")
China3415
#China3415 es un multipolígono:
Out[1]:
| COUNTRY | geometry | |
|---|---|---|
| 0 | China | MULTIPOLYGON (((156893.444 400575.264, 163317.... |
In [2]:
China3415_border = China3415.boundary
type(China3415_border)
Out[2]:
geopandas.geoseries.GeoSeries
In [51]:
China3415_border_gdf = China3415_border.to_frame()
China3415_border_gdf
Out[51]:
| 0 | |
|---|---|
| 0 | MULTILINESTRING ((156893.444 400575.264, 16331... |
In [55]:
China3415_border_gdf.plot(edgecolor="black", linewidth=0.7)
Out[55]:
<Axes: >
Ejercicio 2: Subdivisiones administrativas de China¶
Subdivisiones administrativas de China (1: nivel provincial, 2: nivel de prefecturas)
In [5]:
China_provinces=gpd.read_file(os.path.join("maps","chn_adm_ocha_2020_shp","chn_admbnda_adm1_ocha_2020.shp"))
China_prefectures=gpd.read_file(os.path.join("maps","chn_adm_ocha_2020_shp","chn_admbnda_adm2_ocha_2020.shp"))
Ambas subdivisiones administrativas presentan CRSs no proyectados:
In [6]:
China_provinces.crs, China_prefectures.crs
Out[6]:
(<Geographic 2D CRS: EPSG:4326> Name: WGS 84 Axis Info [ellipsoidal]: - Lat[north]: Geodetic latitude (degree) - Lon[east]: Geodetic longitude (degree) Area of Use: - name: World. - bounds: (-180.0, -90.0, 180.0, 90.0) Datum: World Geodetic System 1984 ensemble - Ellipsoid: WGS 84 - Prime Meridian: Greenwich, <Geographic 2D CRS: EPSG:4326> Name: WGS 84 Axis Info [ellipsoidal]: - Lat[north]: Geodetic latitude (degree) - Lon[east]: Geodetic longitude (degree) Area of Use: - name: World. - bounds: (-180.0, -90.0, 180.0, 90.0) Datum: World Geodetic System 1984 ensemble - Ellipsoid: WGS 84 - Prime Meridian: Greenwich)
Cambiando los CRSs por unos proyectados y que coincidan con los de China3415 para que puedan ser graficados en un mismo mapa:
In [7]:
China_provinces3415 = China_provinces.to_crs(3415)
China_prefectures3415 = China_prefectures.to_crs(3415)
Ejercicio 3: Plantas generadoras de energía en China¶
Datos a emplear: Global Power Plant Database
Es una base de datos que contiene datos geográficos de plantas generadoras de energía de diferentes países. Se creará un GeoDataFrame que incluya solo los datos para China.
In [8]:
import pandas as pd
PowerPlants=pd.read_csv(os.path.join("data","global_power_plant_database.csv"), dtype={'other_fuel3': 'object'})
China_PowerPlants=PowerPlants[PowerPlants["country"]=="CHN"]
China_PowerPlants=China_PowerPlants.loc[:,["name","gppd_idnr","capacity_mw","latitude","longitude","primary_fuel"]]
China_PowerPlants.reset_index(drop=True, inplace=True)
China_PowerPlants.info()
<class 'pandas.core.frame.DataFrame'> RangeIndex: 4235 entries, 0 to 4234 Data columns (total 6 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 name 4235 non-null object 1 gppd_idnr 4235 non-null object 2 capacity_mw 4235 non-null float64 3 latitude 4235 non-null float64 4 longitude 4235 non-null float64 5 primary_fuel 4235 non-null object dtypes: float64(3), object(3) memory usage: 198.6+ KB
In [61]:
China_PowerPlants_gpd=gpd.GeoDataFrame(data=China_PowerPlants.copy(),
geometry=gpd.points_from_xy(x=China_PowerPlants.longitude,
y=China_PowerPlants.latitude),
crs=4326)
#Cambiando los CRSs por unos proyectados y que coincidan con los de China3415:
China3415_PowerPlants=China_PowerPlants_gpd.to_crs(3415)
China3415_PowerPlants.head()
Out[61]:
| name | gppd_idnr | capacity_mw | latitude | longitude | primary_fuel | geometry | |
|---|---|---|---|---|---|---|---|
| 0 | APP Zhenjiang Jindong Mill power station | WRI1075566 | 290.0 | 32.1944 | 119.6998 | Coal | POINT (1047199.039 1756451.972) |
| 1 | Aba Minjiang River Jiangseba | WRI1072501 | 128.0 | 31.4837 | 103.6032 | Hydro | POINT (-502906.088 1699174.246) |
| 2 | Abag Banner Huiteng Liang Phase 1 | WRI1072152 | 49.0 | 43.3500 | 115.9000 | Wind | POINT (667016.352 3043941.695) |
| 3 | Ahai | WRI1000457 | 4750.0 | 27.3488 | 100.5061 | Hydro | POINT (-840063.803 1260296.817) |
| 4 | Aksai A | WKS0070352 | 20.0 | 38.9690 | 94.4210 | Solar | POINT (-1280489.125 2632839.774) |
Se crearán capas para clasificar las plantas acorde a la fuente primaria de energía que emplean:
In [10]:
China3415_PowerPlants["primary_fuel"].value_counts()
Out[10]:
primary_fuel Solar 1318 Hydro 947 Coal 946 Wind 835 Gas 170 Nuclear 12 Oil 5 Geothermal 2 Name: count, dtype: int64
In [50]:
from folium import LayerControl
m = China3415_PowerPlants[China3415_PowerPlants.primary_fuel=='Coal'].explore(color="red",name='Coal',show=False)
m = China3415_PowerPlants[China3415_PowerPlants.primary_fuel=='Gas'].explore(m=m, color="blue",name='Gas',show=False)
m = China3415_PowerPlants[China3415_PowerPlants.primary_fuel=='Geothermal'].explore(m=m, color="green",name='Geothermal',show=False)
m = China3415_PowerPlants[China3415_PowerPlants.primary_fuel=='Hydro'].explore(m=m, color="purple",name='Hydro',show=False)
m = China3415_PowerPlants[China3415_PowerPlants.primary_fuel=='Nuclear'].explore(m=m, color="black",name='Nuclear',show=False)
m = China3415_PowerPlants[China3415_PowerPlants.primary_fuel=='Oil'].explore(m=m, color="orange",name='Oil',show=False)
m = China3415_PowerPlants[China3415_PowerPlants.primary_fuel=='Solar'].explore(m=m, color="gray",name='Solar',show=False)
m = China3415_PowerPlants[China3415_PowerPlants.primary_fuel=='Wind'].explore(m=m, color="brown",name='Wind',show=False)
LayerControl(collapsed=False).add_to(m) #optional
m
Out[50]:
Make this Notebook Trusted to load map: File -> Trust Notebook
Se creará un archivo geopackage que incluya cada capa:
In [863]:
China3415_PowerPlants[China3415_PowerPlants.primary_fuel=='Coal'].to_file(os.path.join("maps","China3415_PowerPlants","China3415_PowerPlants.gpkg"), layer='Coal', driver="GPKG",index="")
China3415_PowerPlants[China3415_PowerPlants.primary_fuel=='Gas'].to_file(os.path.join("maps","China3415_PowerPlants","China3415_PowerPlants.gpkg"), layer='Gas', driver="GPKG",index="")
China3415_PowerPlants[China3415_PowerPlants.primary_fuel=='Geothermal'].to_file(os.path.join("maps","China3415_PowerPlants","China3415_PowerPlants.gpkg"), layer='Geothermal', driver="GPKG",index="")
China3415_PowerPlants[China3415_PowerPlants.primary_fuel=='Hydro'].to_file(os.path.join("maps","China3415_PowerPlants","China3415_PowerPlants.gpkg"), layer='Hydro', driver="GPKG",index="")
China3415_PowerPlants[China3415_PowerPlants.primary_fuel=='Nuclear'].to_file(os.path.join("maps","China3415_PowerPlants","China3415_PowerPlants.gpkg"), layer='Nuclear', driver="GPKG",index="")
China3415_PowerPlants[China3415_PowerPlants.primary_fuel=='Oil'].to_file(os.path.join("maps","China3415_PowerPlants","China3415_PowerPlants.gpkg"), layer='Oil', driver="GPKG",index="")
China3415_PowerPlants[China3415_PowerPlants.primary_fuel=='Solar'].to_file(os.path.join("maps","China3415_PowerPlants","China3415_PowerPlants.gpkg"), layer='Solar', driver="GPKG",index="")
China3415_PowerPlants[China3415_PowerPlants.primary_fuel=='Wind'].to_file(os.path.join("maps","China3415_PowerPlants","China3415_PowerPlants.gpkg"), layer='Wind', driver="GPKG",index="")
In [12]:
China3415_PowerPlants_gpkg = "https://github.com/Ciencia-de-datos-espaciales-2023-2/Evaluacion3_GeoDF-Operaciones/raw/main/maps/China3415_PowerPlants/China3415_PowerPlants.gpkg"
listlayers(China3415_PowerPlants_gpkg)
Out[12]:
['Coal', 'Gas', 'Geothermal', 'Hydro', 'Nuclear', 'Oil', 'Solar', 'Wind']
Ejercicio 4: Provincias al norte y sur de China¶
In [13]:
China3415.centroid
Out[13]:
0 POINT (-451290.166 2334683.574) dtype: geometry
In [14]:
centroidX_CN=China3415.centroid.x[0]
centroidY_CN=China3415.centroid.y[0]
a) Provincias al norte de China¶
In [15]:
China_provinces3415.cx[:,centroidY_CN:].unary_union
Out[15]:
b) Provincias al sur de China¶
In [16]:
China_provinces3415.cx[:,:centroidY_CN].unary_union
Out[16]:
c) Duplicados¶
In [17]:
#Disolver por ADM1_EN (divisiones administrativas a nivel provincial)
China_provinces3415_North = China_provinces3415.cx[:,centroidY_CN:].dissolve(by='ADM1_EN')
China_provinces3415_South = China_provinces3415.cx[:,:centroidY_CN].dissolve(by="ADM1_EN")
In [18]:
# Número de provincias duplicadas:
len(set(China_provinces3415_South.ADM1_PCODE).intersection(set(China_provinces3415_North.ADM1_PCODE)))
Out[18]:
8
In [22]:
base=China_provinces3415_North.plot(facecolor='red', alpha=0.4, edgecolor="white", linewidth=0.3)
China_provinces3415_South.plot(ax=base,facecolor='steelblue', alpha=0.4, edgecolor="white", linewidth=0.3)
Out[22]:
<Axes: >
Ejercicio 5: Convex Hull de plantas de energía solar y nuclear en China¶
a) Plantas de energía solar¶
In [32]:
China3415_PowerPlants_Solar = gpd.read_file(China3415_PowerPlants_gpkg,layer="solar")
#Convex hull
China3415_PowerPlants_SolarHull = China3415_PowerPlants_Solar.dissolve().convex_hull
China3415_PowerPlants_SolarHull
Out[32]:
0 POLYGON ((-26944.862 238706.710, -45202.039 24... dtype: geometry
In [33]:
#Convex hull a GeoDataFrame
China3415_PowerPlants_SolarHull_gdf= gpd.GeoDataFrame(index=[0],
crs=3415,
geometry=China3415_PowerPlants_SolarHull)
China3415_PowerPlants_SolarHull_gdf
Out[33]:
| geometry | |
|---|---|
| 0 | POLYGON ((-26944.862 238706.710, -45202.039 24... |
In [34]:
# Gráfica
base = China3415.plot(facecolor="white", edgecolor="black", linewidth=0.9)
China3415_PowerPlants_Solar.plot(marker='o', color='white',edgecolor='orange', markersize=15,ax=base)
China3415_PowerPlants_SolarHull.plot(ax=base,facecolor='orange',
edgecolor='white',alpha=0.3,
hatch='X')
Out[34]:
<Axes: >
b) Plantas de energía nuclear¶
In [35]:
China3415_PowerPlants_Nuclear = gpd.read_file(China3415_PowerPlants_gpkg,layer="nuclear")
#Convex hull
China3415_PowerPlants_NuclearHull = China3415_PowerPlants_Nuclear.dissolve().convex_hull
China3415_PowerPlants_NuclearHull
Out[35]:
0 POLYGON ((-34893.586 338256.930, -62221.064 58... dtype: geometry
In [30]:
#Convex hull a GeoDataFrame
China3415_PowerPlants_NuclearHull_gdf= gpd.GeoDataFrame(index=[0],
crs=3415,
geometry=China3415_PowerPlants_NuclearHull)
China3415_PowerPlants_NuclearHull_gdf
Out[30]:
| geometry | |
|---|---|
| 0 | POLYGON ((-34893.586 338256.930, -62221.064 58... |
In [36]:
#Gráfica
base = China3415.plot(facecolor="white", edgecolor="black", linewidth=0.9)
China3415_PowerPlants_Nuclear.plot(marker='o', color='white',edgecolor='red', markersize=15,ax=base)
China3415_PowerPlants_NuclearHull.plot(ax=base,facecolor='red',
edgecolor='white',alpha=0.3,
hatch='X')
Out[36]:
<Axes: >
c) Convex hull de plantas de energía solar y nuclear¶
In [37]:
base = China3415.plot(facecolor="white", edgecolor="black", linewidth=0.9)
China3415_PowerPlants_SolarHull.plot(ax=base,facecolor='orange',
edgecolor='white',alpha=0.3,
hatch='X')
China3415_PowerPlants_NuclearHull.plot(ax=base,facecolor='red',
edgecolor='white',alpha=0.3,
hatch='X')
Out[37]:
<Axes: >
Ejercicio 6: Operaciones¶
a) En mapas de China¶
1) Intersección entre mapas de provincias al sur y norte de China = Provincias repetidas
In [38]:
China_provinces3415_SouthNorth_Intersect = China_provinces3415_South.overlay(China_provinces3415_North,
how="intersection", keep_geom_type=True)
China_provinces3415_SouthNorth_Intersect.head()
Out[38]:
| ADM1_ZH_1 | ADM1_PCODE_1 | ADM0_EN_1 | ADM0_ZH_1 | ADM0_PCODE_1 | ADM1_ZH_2 | ADM1_PCODE_2 | ADM0_EN_2 | ADM0_ZH_2 | ADM0_PCODE_2 | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 山东省 | CN037 | China | 中国 | CN | 山东省 | CN037 | China | 中国 | CN | MULTIPOLYGON (((848226.102 2014498.442, 842640... |
| 1 | 陕西省 | CN061 | China | 中国 | CN | 陕西省 | CN061 | China | 中国 | CN | POLYGON ((64694.414 1861474.379, 67245.344 185... |
| 2 | 甘肃省 | CN062 | China | 中国 | CN | 甘肃省 | CN062 | China | 中国 | CN | POLYGON ((-203874.868 1989455.040, -200807.739... |
| 3 | 宁夏回族自治区 | CN064 | China | 中国 | CN | 宁夏回族自治区 | CN064 | China | 中国 | CN | POLYGON ((-203441.498 2261773.633, -212565.869... |
| 4 | 青海省 | CN063 | China | 中国 | CN | 青海省 | CN063 | China | 中国 | CN | POLYGON ((-1541026.129 1943815.359, -1548204.7... |
In [56]:
base = China3415_border_gdf.plot(edgecolor="black", linewidth=0.7)
China_provinces3415_SouthNorth_Intersect.plot(ax=base, facecolor="c", edgecolor="white", linewidth=0.5)
Out[56]:
<Axes: >
2) Diferencia entre el mapa de China y la intersección entre provincias del sur con las del norte
In [40]:
China_provinces3415_Diff = China_provinces3415.overlay(China_provinces3415_SouthNorth_Intersect,
how="difference", keep_geom_type=True)
China_provinces3415_Diff.head()
Out[40]:
| ADM1_EN | ADM1_ZH | ADM1_PCODE | ADM0_EN | ADM0_ZH | ADM0_PCODE | geometry | |
|---|---|---|---|---|---|---|---|
| 0 | Shanghai Municipality | 上海市 | CN031 | China | 中国 | CN | MULTIPOLYGON (((1222926.941 1675247.504, 12249... |
| 1 | Chongqing Municipality | 重庆市 | CN050 | China | 中国 | CN | POLYGON ((72519.004 1698948.257, 75526.165 170... |
| 2 | Zhejiang Province | 浙江省 | CN033 | China | 中国 | CN | MULTIPOLYGON (((1319922.798 1602761.054, 13200... |
| 3 | Jiangxi Province | 江西省 | CN036 | China | 中国 | CN | POLYGON ((753769.308 1119026.001, 755578.523 1... |
| 4 | Yunnan Province | 云南省 | CN053 | China | 中国 | CN | POLYGON ((-369650.354 769812.608, -371355.140 ... |
In [57]:
base = China3415_border_gdf.plot(edgecolor="black", linewidth=0.7)
China_provinces3415_Diff.plot(ax=base, facecolor="c", edgecolor="white", linewidth=0.5)
Out[57]:
<Axes: >
b) En mapas de Brazil¶
Código para los mapas de Brazil del notebook GeoDF_Operations de la clase anterior:
In [42]:
#Brazil, CRS:5641
brazilMaps='https://github.com/Ciencia-de-datos-espaciales-2023-2/geodfprepro/raw/main/maps/brazilMaps_5641_all.gpkg'
brazil=gpd.read_file(brazilMaps,layer='country')
brazil_municipalities=gpd.read_file(brazilMaps,layer='municipalities')
centroidX=brazil.centroid.x[0]
centroidY=brazil.centroid.y[0]
#Municipalidades al sur, norte, sur-norte y centro de Brazil
MunisS_brazil=brazil_municipalities.cx[:,:centroidY]
MunisN_brazil=brazil_municipalities.cx[:,centroidY:]
MunisW_brazil=brazil_municipalities.cx[:centroidX,:]
MunisE_brazil=brazil_municipalities.cx[centroidX:,:]
munisMidWE_brazil=MunisW_brazil.overlay(MunisE_brazil, how="intersection",keep_geom_type=True)
munisMidNS_brazil=MunisN_brazil.overlay(MunisS_brazil, how="intersection",keep_geom_type=True)
muniMidBrazil=munisMidNS_brazil.dissolve().overlay(munisMidWE_brazil.dissolve(), how="union",keep_geom_type=True).dissolve()
#Aeropuertos de Brazil, CRS:5641
Brazil_airports = "https://github.com/Ciencia-de-datos-espaciales-2023-2/Evaluacion3_GeoDF-Operaciones/raw/main/maps/Brazil_airports/airports.gpkg"
airports = gpd.read_file(Brazil_airports)
#Aeropuertos medianos de Brazil
Brazil_AirTopLeft=airports[airports.kind=='medium_airport'].cx[:centroidX,centroidY:]
Brazil_AirTopRight=airports[airports.kind=='medium_airport'].cx[centroidX:,centroidY:]
Brazil_AirBottomLeft=airports[airports.kind=='medium_airport'].cx[:centroidX,:centroidY]
Brazil_AirBottomRight=airports[airports.kind=='medium_airport'].cx[centroidX:,:centroidY]
#Aeropuertos grandes
Large_airport=airports[airports.kind=='large_airport']
D:\anaconda3\envs\Ciencia_De_Datos_Espaciales\Lib\site-packages\geopandas\geodataframe.py:1815: FutureWarning: `unary_union` returned None due to all-None GeoSeries. In future, `unary_union` will return 'GEOMETRYCOLLECTION EMPTY' instead. merged_geom = block.unary_union D:\anaconda3\envs\Ciencia_De_Datos_Espaciales\Lib\site-packages\geopandas\geodataframe.py:1815: FutureWarning: `unary_union` returned None due to all-None GeoSeries. In future, `unary_union` will return 'GEOMETRYCOLLECTION EMPTY' instead. merged_geom = block.unary_union
In [45]:
# Hulls de aeropuertos medianos
Brazil_AirTopLeft_hull=Brazil_AirTopLeft.dissolve().convex_hull
Brazil_AirTopRight_hull=Brazil_AirTopRight.dissolve().convex_hull
Brazil_AirBottomLeft_hull=Brazil_AirBottomLeft.dissolve().convex_hull
Brazil_AirBottomRight_hull=Brazil_AirBottomRight.dissolve().convex_hull
#Hull de aeropuestos grandes
Large_Airport_hull= Large_airport.dissolve().convex_hull
#Gráfica con los convex hulls de aeropuertos medianos
base = brazil.plot(color='white', edgecolor='black')
muniMidBrazil.plot(ax=base,facecolor='magenta',alpha=0.4)
Large_Airport_hull.plot(ax=base,facecolor='black',alpha=0.4)
Brazil_AirTopLeft_hull.plot(ax=base,facecolor='red', alpha=0.4)
Brazil_AirTopRight_hull.plot(ax=base,facecolor='blue', alpha=0.4)
Brazil_AirBottomLeft_hull.plot(ax=base,facecolor='green', alpha=0.4)
Brazil_AirBottomRight_hull.plot(ax=base,facecolor='orange', alpha=0.4)
Out[45]:
<Axes: >
1) Diferencia simétrica entre municipalidades al norte de Brazil con las del centro de Brazil
In [60]:
Brazil_muniMid_SymDiff = MunisN_brazil.overlay(muniMidBrazil, how="symmetric_difference", keep_geom_type=True)
base = brazil.plot(color='white', edgecolor='black')
Brazil_muniMid_SymDiff.plot(ax=base, facecolor="c", edgecolor="white", linewidth=0.2)
Out[60]:
<Axes: >
2) Unión entre municipalidades al sur y este de Brazil (S+E)
In [59]:
Brazil_muniSE = pd.concat([MunisS_brazil,MunisE_brazil],ignore_index=True)
base = brazil.plot(color='white', edgecolor='black')
Brazil_muniSE.plot(ax=base, facecolor="c", edgecolor="white", linewidth=0.2)
Out[59]:
<Axes: >
3) Unión entre aeropuertos medianos al norte de Brazil (TopLeft + TopRight)
In [49]:
Brazil_AirTop = pd.concat([Brazil_AirTopLeft,Brazil_AirTopRight],ignore_index=True)
Brazil_AirTop_hull=Brazil_AirTop.dissolve().convex_hull
#Gráfica con el convex hull resultados de la unión
base = brazil.plot(color="white", edgecolor="black")
Brazil_AirTopLeft_hull.plot(ax=base,facecolor='red', alpha=0.1)
Brazil_AirTopLeft.plot(marker='o', color='white',edgecolor='red',markersize=15,ax=base)
Brazil_AirTopRight_hull.plot(ax=base,facecolor='blue', alpha=0.1)
Brazil_AirTopRight.plot(marker='o', color='white',edgecolor='blue', markersize=15,ax=base)
Brazil_AirTop_hull.plot(ax=base,facecolor="orange", alpha=0.2)
Out[49]:
<Axes: >
